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Abstract 

We investigated the structural and dynamical properties of a tetra- 
hedrally coordinated crystalline ice from first principles based on den- 
sity functional theory within the generalized gradient approximation 
with the projected augmented wave method. First, we report the 
structural behaviour of ice at finite temperatures based on the analysis 
of radial distribution functions obtained by molecular dynamics sim- 
ulations. The results show how the ordering of the hydrogen bonding 
breaks down in the tetrahedral network of ice with entropy increase 
in agreement with the neutron diffraction data. We also calculated 
the phonon spectra of ice in a 3 x 1 x 1 supercell by using the direct 
method. So far, due to the direct method used in this calculation, the 
phonon spectra is obtained without taking into account the effect of 
polarization arising from dipole-dipole interactions of water molecules 
which is expected to yield the splitting of longitudinal and transverse 
optic modes at the T-point. The calculated longitudinal acoustic ve- 
locities from the initial slopes of the acoustic mode is in a reasonable 
agreement with the neutron scatering data. The analysis of the vi- 
brational density of states shows the existence of a boson peak at low 
energy of translational region a characteristic common to amorphous 
systems. 

Keywords: Ice, Density functional theory, Molecular dynamics, Lattice dy- 
namics. 
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1. INTRODUCTION 



Ice, the frozen form of liquid water, is one of the most common materials on 
earth and in outer space, and has important relevance to a large number of 
diver se fields such as astronomy, geo physics, chemical physics, life sciences, 
etc. f Petrenko and Whitworth . Il999| ). Ice-Ih, the hexagonal form of ice is 
the most commonly known phase of ice in which each oxygen atom has four 
neighbours at the corners of a tetrahedron. When water freezes the forces of 
interaction between H 2 molecules win over their thermal motion and they 
form a most stable arrangement precisely with hexagonal symmetry. This is 
the reason why snow flake crystals are always hexagonal. The crystal struc- 
ture of ice at the atomic level marks the way these crystals will look. The 
hydrogen atoms are covalently bonded to the nearest oxygen atoms to form 
water molecules which are linked to each other through hydrogen bonds. 
The high translational symmetry is not retained at the level of the crystal- 
lographic unit cell. The phase diagram of ice is very complex with over 12 
known phases existing, (see Fig. [TJ. Besides the environmental importance, 
ice is also special because of the interesting phenomena contained within its 
structure. The crystal structure of ice is very unusual because, while the 
molecules lie on a regular crystal lattice, there is disorder in their orienta- 
tions. This property leads to many interesting characteristics in electrical 
polarization and conductivity. 

Numerous attempts have been made for the past decade s to understand the 
natur e of the lattice vibrations of ice, in particular, ice-Ih (jMarchi. Tse and Klein , 
Experimental inform ation has been obt ained from infrared absorp tion 



lyoo). Hjxpenmental mrormation nas Deen oDtamea irom mirarea absorption 
( Whallev and Bertiel 1967 ) , Raman scattering ( Klue and Whallev 1 19781) and 



both coherent and incoherent ine lastic neutron scattering (jRenker and Blanckenhagen 
1969; Prask and Trevinol 1972j ). On the theoretical side, three basic ap- 
proaches have been adopted to study the lattice modes. The earliest studies 
involved th e application of lattice dynamics to hypothetical pro ton ordered 
structures ( Whallev and Bertiel . 1967 : Prask and Trevinol . 19721 ). Later lat- 
tice dynam ics was used to study more realis tic orient at ionally disordered 
structures ( Nielsen. Townsend and Ricd . 11984). The most recent work has 



utilized the molecular dynamics simulation techniques (T se. Klein and McDonald , 



19841 ) . Many of the these theoretical studies involved the use of classical mod- 



eled potential t hrough empirical metho d in order to describe the interaction 
of th e system ( Marchi. Tse and Klein . 19861 : Nielson. Townsend and Ricd . 

4J). As result of all these works, the overall features of the lattice mode 
spectrum in the translational region (0-300 cm -1 ) and in the librational re- 
gion (450-950 cm -1 ) are reasonably well understood. Potential-based empir- 
ical modelling has had some success towards the end, but to date, there are 
no empirical potentials capable of describing the ice dynamics and related 
properties across its whole spectra range and describing certain key spectra 
features. 

The ab-initio method has recently gained ground not only because of its relia- 
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bility in the study of static and dyna mics properties of ice (jCote. Morrison. Cui. Jenkins et a! 



ity 

2003r iMorrison and Jenkins! . [1999) but also because the method allo w to 



mode l some important features such as the ordered periodic ice structure flLee. Va nderbi lt. Laarsone 
1993), and the nature of hydrogen bond in different geometries ( Xantheas and Dunning) . 



19931). Our present study to understand the microscopic nature and lattice 



vibratio ns of ice-Ih makes use of the Vienna Ab Initio Simulation Package 
( VASP) (|Kresse and Futhmiiller , 1996J ) which is designed to perform ab-initio 
quantum-mechanical molecular dynamics using pseudopotentials and a plane 
wave basis set and the recent implementation of the projected augmented 
wave (PAW) method. 

The computational study presented in this work begins with the choice of unit 
cell of ice as described in Section 2, which is followed by molecular dynamics 
simulation of ice in a supercell (Section 3) with the unit cell replicated in all 
directions in order to obtain the structural behaviour in comparison with the 
available experimental data. In Section 4 and 4.1 we present the results of 
calculated phonon spectrum of ice in the translational region, which is done 
in a 3 x 1 x 1 supercell in y-direction, as well as the corresponding integrated 
vibrational density of states in Section 4.2. The anomalous behaviour of this 
ice structure observed in the low energy region is presented in Section 4.3. 
The final Summary of this work is presented in Section 5. 



2. COMPUTATIONAL DETAILS 

As mentioned above, the calculations in this work were carried out by using 
the V ienna Ab Initio Simulation Package (VASP) (|Kresse and Futhmiiller 



1996) which has been designed for ab-initio quantum-mechanical molecular 



dynamics simulations using pseudopotentials and a plane wave basis set and 
the recent implementation of projected augmented wave (PAW) method. In 
order to investigate the structural properties of ice, we st arted with the con- 
struc tion of ice crystal using the Bernal and Fowler rule ((Bernal and Fowler , 



19331 ). The rule is based on a statistical m odel of the position of hydrogen 



atoms produced by Pauling ( Pauling . 19351 ) using the six possible configura- 



tions of hydrogen atoms within Ice Ih. It is defined as ideal crystal based on 
the assumptions that: 

• Each oxygen atom is bonded to two hydrogen atoms at a distance of 
0.95 A to form a water molecule; 

• Each molecule is oriented so that its two hydrogen atoms face two, of 
the four, neighbouring oxygen atoms in the tetrahedral coordination; 

• The orientation of adjacent molecules is such that only one hydrogen 
atom lies between each pair of oxygen atoms; 
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• Ice Ih can exist in any of a large number of configurations, each corre- 
sponding to a certain distribution of hydrogen atoms with respect to 
oxygen atoms. 

The schematic drawing shown in Fig. El satisfies one out of the six possible 
orientations of protons of the central water molecule according to these rules. 
Each of the oxygen atoms can be linked to another oxygen by the combination 
of a covalent bond plus a hydrogen bond to form the tetrahedral arrangement 
of oxygen atoms. A unit cell of this ice was prepared in a cubic box according 
to Fig. El with 8 molecules of water. All the atomic degrees of freedom were 
relaxed using VASP with the projected-augmented wave (PAW) formalism at 
high precision. The optimum Monkhorst Pack 4x4x4 fc-point was used in 
addition to the generalized gradient approximation (GGA) of Perdew-Wang 
in order to describe the exchange-correlation and to give good description of 
the hydrogen bonding of water. We used a energy cut-off of 500 eV because 
the 2p valence electrons in oxygen require a large plane wave basis set to 
span the high energy states described by the wavefunction close to the oxygen 
nucleus and also the hydrogen atoms require a larger number of planes waves 
in order to describe localization of their charges in real space. The lattice 
constants of the unit cell were calculated from the plot of the energy against 
the volume. The estimated values of the lattice constants are a = 6.1568 A, 
b = 6.1565 A, c = 6.0816A. The values of a ~ b ^ c imply that the relaxed 
structure is tetragonal with c/a ratio ~ 0.988. The experimental latt ice 
constant reported by Blackman et. al. ( Blackman and Lisgarten . 19581 ) is 



6.35012652 A for the cubic geometry. The final geometry obtained was used 
in our molecular and lattice dynamical study of ice. 



3. MOLECULAR DYNAMICS 

In our molecular dynamics simulation, the final relaxed geometry in Fig. |2] 
(or Fig. El was replicated in all directions using the calculated values of the 
lattice parameters to produce 64 molecules of water which form the hexagonal 
structure shown in Fig. HJ Here the molecular dynamics simulation was done 
for the T-point only in a box of dimension 12.411356 x 12.411356 x 12.259675 
A 3 corresponding to the density 1.01 gem -3 to be compared to the real 
density of ice Ih and ice Ic which is 0.92 gem" 3 . Molecular dynamics was 
carried out for 3 ps using a time step of 0.5 fs. The angle H-O-H of the 
ice structure was compared with liquid water at the different temperatures 
as can be seen in Fig. For ice structure, the H-O-H angle was found to 
be 107.8° compared to liquid water case at different temperature or isolated 
water molecule also calculated with VASP which is found to be 104.5° and 
105.4° degree, respectively, as can be seen in Fig. El The angle shown by solid 
ice is an indication that the oxygen centre preserves its tetrahedral structural 
units. Also, the range of the angular distribution for the high temperature 
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water (or supercritical water) is wider and broader than the corresponding 
liquid water at ambient temperature and the solid ice cases. 
The radial distribution functions for the solid ice accumulated over a 3 ps 
run are plotted at two different temperatures, 100 K and 220 K. The goo 
and goH in Fig. El and ^hh in Fig. Oat 100 K, all exhibit the long-range like 
order when compared to the short range order of liquid water. There is a well 
pronounced first peak of goo a t 2.75 A at 100 K and also the position of the 
first minimum is deeper when compared to the liquid water radial distribution 
functions. This result is co mparable to the result obtained using the classical 
TIP4P modelled potential (jSvishchev and Kusalild . Il994h The result of the 
radial distribution function of oxygen atoms, goo, oxygen- hydrogen atoms 
gon (Fig. [SJ), and hydrogen atoms (Fig. at 220 K were fair ly comparabl e 
with available neutron diffraction scattering data of Soper ( Soper . 2000). 
The position of the first peak (VASP calculation) shows very little difference 
from 100 K while the height of the peak is lower for 220 K due to the effect 
of entropy increase of the hydrogen atoms which results from re-orientation 
of protons that tends to force apart more oxygen atoms to a larger distance. 
The positions of the second minimum for 100 K and 220 K solid ice are, 
respectively, 4.8 A and 4.9 A. The experimental result is slightly shifted to 
the right to the value 5.1 A for 220 K. The deviation from experimental value 
might be due to the phase of crystalline ice under consideration. 



4. LATTICE DYNAMICS 

In this study phonon dispersion cu rves are calculate d by using the PHONON 
package developed by K. Palinski ( Parlinskin . 20021 ) which has been designed 



to take input data of Hellmann-Feymann forces calculated with the help of an 
ab-initio electronic structure simulation program such as VASP. We carried 
out the lattice dynamics study of ice by using the geometry of eight-molecule 
primitive cell discussed in the Section 1. The calculations of force constants 
was carried out by considering a 3 x 1 x 1 supercell containing 24 molecules 
of water which is obtained by matching 3 tetragonal unit cells. At the first 
step of the calculation, the PHONON package is used to define the appropri- 
ate crystal supercell for use of the direct method. As done for the primitive 
unit cell, all the internal coordinates were relaxed until the atomic forces 
were less than 10~ 4 eV/A. The relaxed geometry for a 1 x 1 x 1 supercell 
from the initial configurations containing 8 molecules is shown in Fig. El The 
starting geometry of the molecules in the simulation box shown is such that 
no hydrogen bonds were present but the positions of oxygen atoms follow 
the tetrahedral orientation. After the relaxation, all the protons perfectly 
point to the right direction of oxygen atoms and make the required hydrogen 
bonds necessary as indicated by the dotted lines in Fig. El to preserve the 
tetrahedral orientation of the ice structure. 

Figure HUT I) shows the Brillouin zone belonging to the relaxed structure 
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of our (model) ice shown in Fig. El Figure HUT U) shows the Brillouin zone 
used in the analysis of the measured phonon spectra for the model struc- 
ture of D 2 ice. Let us say once again that the relaxed structure shown in 
Fig. El has the long-range orientational order while the actual structure of ice 
Ih has no long-range orientational order. Therefore, in the analysis of the 
measured phonon dispersion curves of D2O, one uses another model for ice 
shown in Fig. HJ We have to keep this in mind when comparing our calcu- 
lated dispersion curves with the measured frequencies. For the evaluation of 
the phon on dispersion curves we have used t he direct ab-initio force constant 
method ( Parlinskin. Li and Kawazoel . 1997 ). whereby the forces are calcu- 



lated with VASP via the Hellmann-Feymann theorem in the total energy 
calculations. Usually, the calculations are done for a supercell with periodic 
boundary conditions. In such a supercell, a displacement u(0, k) of a single 
atom induces forces F(lk) acting on all other atoms, 

F a (lk) = J2® al 3(lk;l f k')u ( i(l'k>). (1) 

Vk'P 

This expression allows to determine the force c onstant matrix directly from 
the calculated forces ( see Parlinski et. al.) ( Parlinskin. Li and Kawazoel . 



19971 IParlinskml 2002). The phonon dispersion branches calculated by the 



direct method are exact for discrete wave vectors defined by the equation 

exp (2mk L • L) = 1, (2) 

where L = (L a , Lb, L c ) are lattice parameters of the supercell. A related tech- 
nique has recently been used to obtain accurate full phonon dispersions in 



highl y symmetric structures of Ni 2 GaMn (jZavak. Entel. Enkovaara. Avuela et al 
2003J. 



In order to obtain the complete information of the values of all force con- 
stants, every atom of the primitive unit cells was displaced by 0.02 A in both 
positive and negative non-coplanar, x, y and z directions to obtain pure har- 
monicity of the system. We use a 3 x 1 x 1 supercell which implies that 
3 points in the direction [100] are treated exactly according to the direct 
method. The points are [C00], with ( = 1, 1/3, 2/3. We calculate forces 
induced on all atoms of the supercell when a single atom is displaced from 
its equilibrium position, to obtain the force constant matrix, and hence the 
dynamical matrix. This is then followed by diagonalization of the dynam- 
ical matrix which leads to a set of eigenvalues for the phonon frequencies 
and the corresponding normal-mode eigenvectors. The vibrational density of 
states (VDOS) is obtained by integrating over fc-dependent phonon frequen- 
cies from the force-constant matrix in supercells derived from the primitive 
molecule unit cells. 
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4.1 Phonon dispersions of crystalline ice 



The phonon dispersion curves calculated for our ice crystal in [100] direc- 
tion are shown in Fig. ^JJ for low lying energy vibrations. According to 
the geometry of the supercell, the low-frequency (or low-energy 0-50 meV) 
acoustic mode s can b e com pared to Renker's inelastic neutron scattering 
measurement ( Renkerl . Il973l ) along the [0001] direct ion (T A) of hexagonal 
symmetry shown in Fig. EH taken from reference ([hi Il996h . though our cal- 
culation was done only along the Cartesian direction [100] of the cubic cell. 
Our transverse and longitudinal acoustic (LA/TA) dispersions are well be- 
haved when compared to some other modelled calculations or experimental 
results ( Renkerl . Il973l ) in the high symmetry directions TA of hexagonal ice 
as shown in Fig. ITTTc). We can also compare our result to the Cote et. al. 
result in Fig. ITTT b) (|Cote. Morrison. Cui. Jenkins et ~al. . 2003f ). where they 
have recently used the ab-initio method to obtain the phonon dispersions in 
the translational frequency range for the ice structure in the Brillouin zone 
of the orthorhombic eight-molecule unit cell. Altogether our LA and TA dis- 
persions are better than Cote's LA/TA in comparison to the experimental 
curves in Fig. ITTT c) . Our dispersion curves in [100] direction can as well be 
compared to the dispersion curves obtained using a dynamical model with 
two force constants to des cribe the low frequencies of vibrations of hexagonal 
ice as proposed by Faure ( Faurel 1969I ). 



Although everywhere along the T-point, our dispersions are completely de- 
generate in the optic region whereas Cote's dispersions in Fig. ITTT b) show 
some splittings, the so-called longitudinal/transversal optic (LO/TO) split- 
tings whose origins is explained below, while at the zone boundary, some of 
the dispersions are non-degenerate unlike our results. Our inability to re- 
produce these splittings at the T-point is due to the direct method approach 
in which absolute periodicity of the crystal according to Born-von Karman 
conditions was considered. The splitting of LO and TO branches for long 
wavelengths occurs in almost all crystals which are heteropolar (partially 
ionic such as GaAs) or ionic (such as NaCl) at the T-point, and only for in- 
frared active modes (|Srivastava . 1990). The long-range part of the Coulomb 
interaction causes the splitting of the k = optic modes raising the frequency 
of LO modes above those of TO modes. The long-range part of the Coulomb 
interaction corresponds to the macroscopic electric field arising from ionic 
displacements. Ice is a tetrahedrally covalently bonded polar system whose 
dipole-dipole interactions give rise to the electric field when they are dis- 
turbed. The origin of the splitting is therefore the electrostatic field created 
by long wavelength modes of vibrations in such crystals. Usually a micro- 
scopic electric field influences only the LO modes while TO modes remain 
unaltered. The field therefore breaks the Born-von Karman conditions, as a 
consequence with a direct method only finite wave vector k ^ calculations 
are possible. Elongated su b-supercells are n eeded to recover the k — > limit 
of the LO phonon branch ( Parlinskin . 20021 ). 
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In our result there are two transverse acoustic branches which are highly 
degenerate and a longitudinal acoustic branch. The first optical branch of 
the dispersion curves is degenerate with the transverse acoustic branches 
at energy ~9.0 meV. The transverse and longitudinal velocities of sound are 
calculated from the initial slopes of the corresponding transverse and longitu- 
dinal acoustic branches of the dispersion curves in the long wavelength limit. 
The experimental values of velocities reported in Tableware those of longitu- 
dinal and transverse sound waves propagating along the c-direction of single 
crystals of ice at 257 K. It is well known that velocities of sound depend much 
on the direction of propagation and also on the temperature. Inelastic X-ray 
scattering data from water at 5 °C shows a variation of the sound velocity 
from 2000 to 3200 m/s in the momentum range of 1-4 nm' 1 . The so-called 
transition from normal to fast sound in liquid water at ~ 4 meV, the energy 
of sound excitations which is equal to the observed second weakly dispersed 
mode, was reported to be due t o the reminiscent of a phonon branch of ice 



Ih of known optical character (jSette. Ruocco. Krisch. Masciovecchio et al. 



1996). We can conclude that our calculated values of longitudinal velocity, 
vl is in a reasonable range of velocity of sound in ice along the [100]-direction 
chosen for our calculation. We must also stress the fact that our phonon dis- 
persions were calculated at K. 



4.2 Vibrational density of states of crystalline 
ice 

In order to understand the mode of collective vibration of molecules of wa- 
ter in ice from a spectroscopic point of view, we need to consider the three 
normal modes of an isolated water molecule shown in Fig. as phases of 
water vapour; liquid water and ice consist of distinct H 2 molecules rec- 
ognized by Bernal and Fowler in 1933, which explained one of the factors 
leading to the Pauli model of the crystal structure of ice. The fact that 
the forces between the molecules are weak in comparison with the inter- 
nal bonding results in a simple division of the lattice modes into three 
groups involving the internal vibrations, rotations, and translations of the 
molecules. The frequency of the first two groups depend primarily on the 
mass of the hydrogen or deuterium nuclei, and th e frequencies of the trans- 



lation s depend on the mass of the whole molecule ([Petrenko and Whitwo rth. 



1999). A free H2O molecule has just three normal modes of vibration il- 



lustrated in Fig. The comparatively small motions of the oxygen atoms 
are required to keep the centre of mass stationary, and these motions re- 
sult in the frequency u 3 being slightly higher than v\\ these depend on 
the force constant for stretching the covalent O-H bond, while the bend- 
ing mode v-i depends on the force constant for changing the bond angle. In 
the vapour the free molecules have a rich rotation-vibration infrared spec- 



S 



trum (jBenedict. Gailar and Plvlerl . 1195 6). from which the frequencies of the 



molecular modes are deduced to be: 



v\ = 3656.65 cm 1 = 453.4 meV, 
i/ 2 = 1594.59 cm" 1 = 197.7 meV, 
zy 3 = 3755.79 cm" 1 = 465.7 meV. 

For ice the band around 400 meV is thus ratified with the O-H bond stretch- 
ing modes V\ and v^. The frequencies are thus lowered from those of the free 
molecules by the hydrogen bonding to the neighbouring molecules, but as 
a single molecule cannot vibrate independently, this coupling also leads to 
complex mode structures involving many molecules. 

We can now discuss the vibrational density of states, VDOS, for H 2 ice 
based on the lattice dynamics obtained from the results of our calculation 
and compare them to some of the well known spectra of ice such as in- 
frared and Rahman spectra and inelastic neutron scattering data. The total 
VDOS calculated from the phonon dispersions for our ice structure is shown 
in Fig. Also shown in Fig. ED is the corresponding partial VDOS for 
both hydrogen and the oxygen atoms in the ice system. We note that these 
phonon DOS are not complete since the summation is not done over the 
whole Brillouin zone, but only in the [100] direction of the cubic symmetry. 
The distribution of the partial DOS is given by 

ffa,*H = ^ M^; k >i)l 2( W^ -^(k,j)), (3) 

where e a (k; k, j) is the a-th Cartesian component of the polarization vector 
for the fc-th atom; n is the number of sampl ing points and d is the dimension 
of the dynamical matrix ( Parlins kinl . l2002h . The total VDOS is calculated 



by summing all the partial contributions. Figure [TBI shows the total VDOS 
together with the full phonon dispersion curves along the [100] direction. 
Also shown in Fig. [T3| is the enlargement of the intermolecular frequency 
range on which we su perimpose the inelastic neutron-scat t ering spectra data 



extracted from Ref. ( Cote. Morrison. Cui. Jenkins et al . |2003). The com- 



parison is made with the ice Ih data though our ice structure is not perfectly 
hexagonal but still there is very little difference between the neutron data 
for ice Ih and ice geometry used in our calculation in the translational region 
as explained below. There are well defined separated peaks in the whole 
range of the vibrations. The illustrative discussion in Fig. ^] can be well 
understood if we consider the partial DOS in Fig. The covalently O-H 
stretching mode of both phase and anti-phase, analogous to the frequencies 
v\ and ^3 for isolated free water molecules, can be seen clearly in Fig. IToT b) 
in the energy range (350-410 meV) or frequency range (3010-3400 cm -1 ). We 
can notice that the collective motion of oxygen is almost static when com- 
pared to the collective contributions from the hydrogen atoms. According to 
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the Ra hman spectra, a strong peak is observed at 382 . 3 meV (3083 cm 1 at 
95 K) ([Bertie and Whallevl . fl967t IWhallev and Bertiel . Il967i ) for D 2 0. If we 
take into account the mass difference between deuterium and hydrogen atoms 
(i.e., isotope effect), the peaks which are observed at 2950, 3000 (very short), 
3250, and 3270 cm -1 are in good range when compared to the experimentally 
observed values for v\ and z/3. In the intra-molecular bending region, anal- 
ogous to the frequency u 2 for an isolated water molecule (1580-1680 cm -1 ), 
there is an interesting feature. Our results show that, of all the contributions 
resulting from the collective motion of hydrogen atoms as contribution from 
collective motion of oxygen atoms is recessive, only one of the components of 
the collective motions of hydrogen atoms contributes to the intra-molecular 
bending modes and it is one that is dominant. Figure flBT b) gives an ex- 
ample of such contribution being dominated mainly by the ^/-component of 
the intra-molecular vibration of the O-H. This means that intra-molecular 
bending of the angular motion takes place mostly in one direction. 
Tanaka has identified hydrogen-bond bendi ng modes with negative expan- 
sion coefficients associated with this region (jTanakal Il998h . If we go further 
down to the low frequency region such as 600-1200 cm -1 called the molecu- 
lar librational region and then to (0-400 cm -1 ) called the molecular transla- 
tional region, where we estimated the sound velocities from the corresponding 
phonon dispersion curve, the VDOS peaks of these modes of vibration agree 
very well with the experimental observation from inelastic neutron scattering 
data. The general agreement of the features in the translational optic region 
i s good w ith all the three distinct peaks present at 400, 270 and 105 cm -1 
(p. ll99fih . 



4.3 The boson peak in ice 

The high frequency (0.1-10 THz) or energy (0.4136-41.36 meV) excitations 
have been experimentally shown to have linear dispersion relations in the 
mesoscopic momentum region (~ 1-10 nm _1 ). So many amorphous materi- 
als display low temperature anomalies in their specific heats that they are 
generally regarded as being universal properties of the glassy state. These 
anomalies are usually of two kinds. The first concerns observation that, while 
many crystals obey the Debye law C oc T 3 for temperatures less than say, 
IK, glasses with the same chemis try frequently display the law C oc T at 
correspondingly low temperature (|Phillipsl . 1996). This second observation 



is tied in with the appearance of the ubiquitous Boson peak (BP) in inelasti c 



neutron and Raman spectra ([Buchenau. Zhou. Nucker. Gilrov et all 119871 ) . 



The name Boson peak therefore refers to the fact that the temperature de- 
pendence of its intensity scales roughly with the Bose-Einstein distribution. 
Dove et. al. demonstrated that the Boson peak arises largely from a flatten- 
ing of the dispersion of the transverse acoustic modes, in a manner similar 
to that which occurs in a crystalline material at its Brillouin zone bound- 
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ary (jDove. Harris. Harmon. Parker et al 



(Elliott, 199 



/ suggested by | 
19961 iTaraskint 



original l y suggested bv fl Leadbcttcr. 1969) , and recently de veloped further by 



19971 ). This idea is not new, being 



19971 ). A standard way 



Il997t iTaraskin and Elliot , 
of extracting the boson peak (BP) from the vibrational spectrum is to plot 
g(E) j E 2 (as done in Fig. llfif b)). since in the Debye approximation g(E) ~ E 2 
at low energy. In Fig. EI] we show the plot of g(E) and the corresponding 
g(E) j E 2 in the translational low energy range for the ice geometry in our cal- 
culation. According to the Debye law, it is expected that g(E)/E 2 should be 
constant for the whole range of energy. This constant relation is only obtained 
at energies larger than 40 meV in agreement with the experimental obser- 
vation range for the BP as mentioned above. There is an anomalous sharp 
peak at 3.5 meV which can be ascribed to the region of low-energy excess vi- 
brational excitation of the so-called BP. The reason for this peak is unknown 
since we have a crystalline structure for our ice geometry. The peak reveals 
the anomalous behaviour of hydrogen bonding in the crystal ice which shows 
similarity in the behaviour as in the results of an inelastic neutron scattering 
study of a crystalline polymorph of SiC>2 (a-quartz), and a number of silicate 



glasses (pure silica, Si0 2 ) with tetrahedral coordination ( Harris et al . 2000[ ). 



Also amorphous solids, mo st supercool liquids and the complex sys t ems s how 
this anomalous character ( Grigera. Martin. Parisi and Verrocchiol . 2003). 



5. SUMMARY 

Our MD simulation, carried out through the analysis of radial distribution 
functions of the ice crystal, shows fair agreement in the positions of peaks in 
comparison with neutron diffraction data. 

The phonon dispersion calculations in [100] direction shows a good result 
especially in the acoustic region in comparison with experiment. Our optic 
modes are degenerate at the T-point due to application of the direct method 
used in this calculation which does not take into account the effect of po- 
larization arising from dipole-dipole interactions of water molecules which 
is expected to yield a significant effect in the splitting of longitudinal and 
transverse optic modes at the T-point. We intend to include this polariza- 
tion effect in our future work through the calculation of effective charge ten- 
sors using the Berry phase approach to see actually if there will be splitting. 
Our calculated longitudinal acoustic velocity agrees well with the longitudi- 
nal acoustic velocity from inelastic neutron scattering data. The vibrational 
density of states reproduces all the features in covalently O-H stretching re- 
gion, intra-molecular bending region, molecular librational region as well as 
in the molecular translational region, when comparing our results to some 
infrared spectra, Rahman spectra as well as inelastic neutron scattering re- 
sults. 

The analysis of the vibrational density of states shows a boson peak, a char- 
acteristic feature common to amorphous systems, at low energy of the trans- 
lational region. The anomalous sharp peak of g(E)/E 2 at 3.5 meV, which 
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can be ascribed to the region of low-energy excess vibrational excitation of 
the boson peak, might be due to the anomalous behaviour of hydrogen bond- 
ing since we have a perfect crystal. The origin of this anomalous behaviour 
of this abnormal peak in non-crystalline solids and supercool liquids is still 
subject of scientific debate. 
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TABLES 

Table U Velocities of sound calculated from the initial slope of the phonon 
dispersion curves of ice in [100] direction compared to the exper- 
imental result. 
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FIGURE CAPTIONS 



Fig. ^ The phase diagram of the stable phases of ice. 



Fig. 121 
Fig. El 



Fig. H 



Fig.[U 



Fig. El 



Fig.m 



Fig. Ill 



The tetrahedral unit from which the hexagonal ice is created. 

Initial and the relaxed geometry of the unit cell of ice. The ice 

structure was initially packed in a cub ic unit cell with initial lat- 

tice constant taken from the literature (Lee. Vanderbilt. Laarsonen. Car et al 



1993) to be 6.35 A. There are no hydrogen bonds in the initial pre- 



pared structure shown on the left but were perfectly formed after 
the relaxation. The relaxed geometry has the values of a pa b ^ c 
which implies that the relaxed structure is tetragonal with c/a 
ratio pa 0.988. 

Relaxed crystalline structure of ice produced from the geometry 
in Fig. |21 by replicating in all direction by the calculated lattice 
parameters. 

Distribution of H-O-H angle in ice and water (in arbitrary units). 
The comparison is done for the ice at 220 K and liquid water 
simulated at room temperature (298 K) and at high temperature 
in the supercritical regime. The H-O-H angles of ice structure are 
larger compared to those of liquid water. 

Radial distribution functions goo and gon at 100 K and 220 K 
obtained with VASP. There is a gradual loss of long-range order 
at 220 K as can be noticed in its 2nd and 3rd peaks of goo when 
compared to the results for 100 K. 

Radial distribution functions ^hr" at 100 K and 200 K obtained 
with VASP. There is a gradual "loss of peaks" observed for tem- 
perature 100 K at 2.5, 4.4, 4.8 and 5.4 A due to the loss of long- 
range order as the temperature increases. 



Radial distribution functions goo and goa at 220 K obtained with 
VASP compa red to neutron diffraction scattering data (NDS) 
(jSoperL l200m . 



Fig. EH Radial distribution functions #hh obtained with VASP at 220 K 
comp ared to neutron diffraction scattering data (NDS) (|Soper . 
l2000h . 

Fig. 1101 (I) Brillouin zone for the cubic symmetry used in our VASP cal- 
culation of the ice system. Our phonon dispersions are calcu- 
lated in [100] /^-direction, (II) The first Brillouin zone for the 
structure of ice Ih with origin at the point T. TA =\c* and TM 
=la*, where a* and c* are the vectors of the reciprocal lattice 
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( Pc trenko and Whitworthl . 1999). The dispersion curves are com- 
monly drawn along the lines of symmetry TA, TM, and TK. 

Fig. 1111 (a) With VASP calculated phonon dispersion curves of ice in [100] 
direction of the tetragonal unit cell compared to (b) dispersion 
relations in the translational frequency range for the ice structure 
plotted along the Cartesian directions from zone center to zone 
edge in the Brillo uin zone of the orthorhombic eight-molecul e unit 
cell (Cote et. al. ( Cote. Morrison. Cui. Jenkins et al. . 2003h ) and 
(c) the experimental d ispersion of D2O ice according to Renker's 
model ((Renker . 1973f ). The difference in scale of (c) from (a) and 
(b) is due to the isotopic effect because of difference in the masses 
of hydrogen and Deuterium atoms. 



Fig. UJ 



Fig.UJ 



Fig. m 



Fig. [[5] 



Fig. [U 



The three normal modes of an isolated water molecule. Mo- 
tion with frequency v\ can be regarde d as symmetric stretch, vo 
as be nding and z/3 as anti-symmetric (jPetrenko and Whitworth . 
IT" 



Total vibrational density of states (VDOS) of ice based on the 
lattice dynamics together with full phonon dispersion curves. The 
VDOS show all the important regions such as the intermolecular 
translational, librational, bending and the stretching frequency 
range. 

Enlargement of calculated total vibrational density of states in 
Fig. II 31 showing the intermolecular range of frequencies. The bro- 
ken line is taken fr om the inelastic neutron scattering data a vail- 



able through Ref. (jCote. Morrison. Cui. Jenkins et al. , 2003). 



Partial density of states of ice based on the lattice dynamics, (a) 
shows the x, y and z components of VDOS for the oxygen atoms 
in the H 2 ice. (b) shows the x, y and z VDOS for corresponding 
hydrogen atoms for the whole range of frequencies. It is interest- 
ing to notice that in the intra molecular bending region (1580-1680 
cm -1 ), only one of the components of the VDOS of hydrogen (say 
y) dominates while x and z contribute less. The contribution from 
the oxygen atoms can be regarded as being completely recessive 
in this region. 

Plot of the VDOS g(E) and the corresponding g(E)/E 2 vs. E for 
the region of translational mode. The boson peak is found in the 
low-energy region at 3.5 meV. 
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Table 1: W. A. Adeagbo et. al. 
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a (Gammon, Kicftc. Cloutcr and Dcnner, 1983) 
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In each tetrahedral unit there 
are 6 different ways to arrange 
hydrogen 



Figure 2: W. A. Adeagbo et. al. 
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Figure 8: W. A. Adeagbo et. al. 
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Figure 11: W. A. Adeagbo et. al 
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Figure 12: W. A. Adeagbo et. al. 
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Figure 13: W. A. Adeagbo et. al. 
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Figure 14: W. A. Adeagbo et. al. 
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Figure 15: W. A. Adeagbo et. al. 



25 




10 20 30 40 50 60 
E (meV) 



Figure 16: W. A. Adeagbo et. al. 



26 



